setwd("/home1/30/jc227089/SRE_inv/ran_outs")

source("/home1/30/jc227089/evo-dispersal/SRE_inv/SRE_inv_functions.R")

#mother function for running races - keeps space same length and returns time to reach end of space
# can do spread up to 50
mother<-function(n=20*40*10, ngens.init=100, gens.breach=1000, spY=20, steps, K=40, lambda, d.cost, mutn=0.05, plot=F, dval=99, spread=0){
	pop<-init.inds(n, 10, spY, hvar=T, dval)
	if (dval<=1) dmut<-FALSE
	else dmut<-TRUE
	maxX<-75
	hab<-init.hab(10, spY, 0, spread=0)	
	nbrs<-neighbours.init(10, spY) #reveal only x up to 10
	for (i in 1:ngens.init){
		pop<-repro.disp(pop, hab, 10, spY, K, lambda, d.cost, nbrs, mutn, adap=T, dmut=dmut)
		if (plot==T) plotter(pop, maxX, spY, K)
	}
	init.traits<-col.sample(pop, 5) #sample initial means and variances
	nbrs<-neighbours.init(maxX, spY) #reveal all neighbours
	ngens<-ngens.init+spread+gens.breach
	hab<-init.hab(maxX, spY, steps, spread=spread)
	st.time<-0
	for (i in (ngens.init+1):ngens){
		pop<-repro.disp(pop, hab, maxX, spY, K, lambda, d.cost, nbrs, mutn, adap=T, dmut=dmut)
		if (plot==T) plotter(pop, maxX, spY, K)
		if (max(pop[,"X"])==maxX | i==ngens) {
			return(i)
		}
	}
}

# initialises a habitat matrix based on homogenous area (maxX-nsteps, spY) and nsteps
# nsteps defines the rate at which habitat declines beyond homogenous area and >0
# homogenous area has default survival value of 1
# gradient starts 10+spread units from LHS
init.hab<-function(maxX, spY, nsteps, hab.val=1, spread){
	if (nsteps>maxX) return(print("Error, maxX>nsteps"))
	if (nsteps==0) return(matrix(data=1, nrow=maxX, ncol=spY))
	ncol1<-10+spread
	ncol2<-maxX-10-spread-nsteps
	simple<-matrix(hab.val, nrow=spY, ncol=ncol1)
	simple2<-matrix(0, nrow=spY, ncol=ncol2)
	d<-seq(hab.val, 0, length.out=nsteps+1)[-1]
	d<-rep(d, times=spY)
	grade<-t(matrix(d, ncol=spY, nrow=nsteps))
	t(cbind(simple, grade, simple2)) #same format as density matrix
}


### initialise ###
args=(commandArgs(TRUE))

#evaluate the arguments
# input arguments will be number of replicates runs (rr), file.ID, spread time (sp)
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}


sim.result<-c()
for (rep in 1:rr){
	out<-mother(steps=10, lambda=5, d.cost=0.1, plot=F, spread=sp, gens.breach=3000)
	sim.result<-rbind(sim.result, out)	
}
sim.result<-cbind(rep(sp, rr), sim.result)
save(sim.result, file=paste("SRE_inv_race", file.ID, ".RData", sep=""))
